JWST/MIRI Simulated Imaging: Insights into Obscured Star Formation and AGNs for Distant Galaxies in Deep Surveys

您所在的位置:网站首页 false `constrain JWST/MIRI Simulated Imaging: Insights into Obscured Star Formation and AGNs for Distant Galaxies in Deep Surveys

JWST/MIRI Simulated Imaging: Insights into Obscured Star Formation and AGNs for Distant Galaxies in Deep Surveys

2023-03-22 02:17| 来源: 网络整理| 查看: 265

1. Introduction

Mid-infrared (mid-IR) wavelengths provide extremely valuable information for extragalactic sources. Galaxies often contain a large amount of polycyclic aromatic hydrocarbon (PAH) molecules (see Tielens 2008 for a review). When heated by starlight, these PAH molecules produce strong emission features mainly at ≈3–18 μm. In star-forming galaxies, PAHs can emit as much as ∼20% of the total IR luminosity, and the 7.7 μm PAH feature can contribute as much as 50% of the total PAH emission (e.g., Shipley et al. 2016). The PAH emission features not only can reveal important information about, e.g., metallicity and ionization parameter (e.g., Draine & Li 2007; Shivaei et al. 2017) but also are potential redshift indicators for high-redshift objects (e.g., Chary et al. 2007).

Active galactic nuclei (AGNs) are often surrounded by large amounts of dust (e.g., Antonucci 1993; Urry & Padovani 1995; Netzer 2015). The dust absorbs UV/optical radiation from the central engine and reaches temperatures of above a few hundred kelvins, much higher than the typical temperatures of interstellar dust heated by starlight (a few tens of kelvins). This AGN-heated hot dust reemits the absorbed energy mainly at mid-IR wavelengths. Therefore, AGNs can be identified based on the hot-dust emission (e.g., Stern et al. 2012; Assef et al. 2013; Kirkpatrick et al. 2015). The mid-IR selection of AGNs has significant advantages over other methods in terms of dust obscuration (e.g., Hickox & Alexander 2018; Alberts et al. 2020). The UV/optical selection is often affected by dust obscuration and biased to bright type 1 AGNs. X-ray data are currently the most robust method of AGN selection (e.g., Brandt & Alexander 2015; Xue 2017). Deep X-ray surveys can sample low-luminosity AGNs below LX ∼ 1043 erg s−1 up to z ∼ 4 (Xue et al. 2016; Luo et al. 2017; Vito et al. 2018). At such a low LX level, most (≈90%) of the X-ray-selected AGNs are type 2, which are potentially missed by UV/optical selections (e.g., Merloni et al. 2014). However, X-ray selection could miss a large population of extremely obscured "Compton-thick" AGNs (neutral hydrogen column density, NH > 1024 cm−2; e.g., Brandt & Alexander 2015; Hickox & Alexander 2018). These obscured AGNs can be identified through warm-dust emission from the AGNs in the mid-IR (e.g., Alexander et al. 2008; Del Moro et al. 2016), and therefore it is useful to explore the utility of mid-IR observations to identify these objects.

Past and current facilities have not been sufficiently sensitive to capture the mid-IR fluxes from most sources in the distant universe. Until now, the most sensitive mid-IR facility has been Spitzer/IRAC+MIPS, covering wavelengths 3.6–24 μm. However, the IRAC+MIPS coverage has a large uncovered gap between the longest IRAC band (8 μm) and the shortest MIPS band (24 μm), leaving a lot of room for model degeneracy at mid-IR wavelengths. For example, observations have found some z ∼ 2 galaxies with elevated 24 μm emission compared to that expected from star formation (mid-IR-excess galaxies; e.g., Daddi et al. 2007b; Papovich et al. 2007). The mid-IR excess emission may be interpreted as either PAH emission or AGN-heated dust radiation, and it is challenging to decompose these two components with broadband imaging alone that lacks contiguous mid-IR wavelength coverage (e.g., Daddi et al. 2007a; Azadi et al. 2018). Furthermore, the Spitzer imaging has a large point-spread function (PSF; FWHM ≈ 6'' at 24 μm), causing source confusion issues in crowded deep fields (e.g., Dole et al. 2004; Ashby et al. 2018). AKARI/IRC provided contiguous wavelength imaging over 2–26 μm. However, its sensitivity is relatively low compared to Spitzer, and many Spitzer faint sources are undetectable by AKARI (e.g., Papovich et al. 2004; Clements et al. 2011). Due to the lack of sensitive 8–24 μm imaging, many of the IR AGN selection techniques are forced to use ≲8 μm colors that are also biased toward type 1 AGNs, similar to UV/optical selections (e.g., Donley et al. 2012; Li et al. 2020).

The upcoming James Webb Space Telescope (JWST) mission will revolutionize the field of infrared astronomy. The onboard Mid-Infrared Instrument (MIRI), equipped with both an imaging camera and a spectrograph, covers a wavelength range from ≈5 to 28 μm. The sensitivity of MIRI is generally ≈1–2 orders of magnitude higher than Spitzer (e.g., Glasse et al. 2015), able to capture the bulk of cosmic mid-IR emission (e.g., Bonato et al. 2017; Rieke et al. 2019). The MIRI imager has a PSF FWHM of ≃0farcs2–0farcs8. The subarcsecond angular resolution provided by MIRI reduces greater source blending and source confusion. The detector's field of view (FOV) is $1\buildrel{\,\prime}\over{.} 3\times 1\buildrel{\,\prime}\over{.} 7$, sufficient for extragalactic surveys. MIRI has a total of nine broadband filters with contiguous wavelength coverage and is able to characterize well the mid-IR spectral shape.

Aside from observational facilities, reliable techniques are also essential for the estimation of source properties. Past methods have relied on color selection (and color–color selection) to classify IR sources (e.g., Stern et al. 2005; Donley et al. 2012; Messias et al. 2012; Kirkpatrick et al. 2017). The color–color methods provide a mechanism to classify galaxies as AGN versus non-AGN, and to identify candidates for "composite sources" (e.g,. Kirkpatrick et al. 2017). These methods can be be straightforwardly applied to a large sample of sources. Indeed, previous work has shown that MIRI mid-IR color–color diagrams can reliably identify AGNs in galaxies down to low AGN luminosities (Eddington ratios of ∼0.01; e.g., Kirkpatrick et al. 2017). However, color–color methods have disadvantages. They often require knowledge of the redshift of the source a priori to interpret (e.g., the redshift can determine which colors best separate warm dust heated by AGNs from emission from PAHs). In addition, they are typically binary (either a source is an AGN or it is not), or trinary (including composite sources). Lastly, typical color selection methods use only a portion of the available information (e.g., a color–color diagram only uses flux densities measured at ≤4 passbands). Historically, these methods have been very successful. However, as the amount of available data and passbands increases (as will be the case in the JWST/MIRI era), it is prudent to explore methods that can provide more detailed constraints using the full data set and thereby improve our understanding of the coevolution of AGNs and star formation in distant galaxies.

An alternative type of quantitative method to characterize the IR emission from sources is spectral energy distribution (SED) fitting. SED fitting has been successfully applied to AGN identification/characterization in previous works, based on the fact that AGNs and star-forming galaxies typically have distinctive SED shapes (e.g., Alonso-Herrero et al. 2006; Caputi 2013; Chang et al. 2017). SED-fitting methods typically fit the photometric data with numerous model templates. They are more computationally intensive than color–color methods, but they have the advantage that they can simultaneously take into account all photometric data and fit for multiple model parameters, including fitting for a photometric redshift (e.g., Małek et al. 2014). Simulations indicate that MIRI data, when included in the SED fitting, can significantly improve the accuracies of fitting results (e.g., Bisigello et al. 2016, 2017; Kauffmann et al. 2020).

Here we employ one such SED-fitting method, x-cigale, which is an efficient python code that can model multiwavelength galaxy/AGN SEDs from X-ray to radio (Yang et al. 2020). Compared to its original version, cigale (Boquien et al. 2019), x-cigale implements many AGN-related improvements. For example, x-cigale allows a polar-dust (PD) component that has been recently observed with high-resolution imaging of local AGNs (e.g., Asmus 2019; Stalevski et al. 2019). Thanks to its optimized parallel algorithm, x-cigale is able to fit thousands of SEDs using millions of models within a few hours on a typical multicore desktop/laptop. x-cigale adopts physical models following the law of energy conservation. Aside from the traditional least-χ2 analysis, x-cigale also allows simultaneous analysis for source properties such as redshift and fracAGN (where fracAGN is the fraction of total dust IR luminosity attributed to the AGN) in a Bayesian style. The Bayesian analysis considers the full probability density functions (pdf's). In contrast, the least-χ2 analysis only considers the best-fit SED model, and the nonnegligible probability of other models may be neglected. Therefore, the Bayesian results provide estimates of the pdf's for model parameters that are more informative than minimum-χ2 results (e.g., Pirzkal et al. 2012; Boquien et al. 2019).

In this paper, we investigate the potential ability of deep imaging with JWST/MIRI to constrain the properties of distant galaxies. We use as our fiducial example the MIRI observations expected as part of the Cosmic Evolution Early Release Science (CEERS) survey, including the existing multiwavelength photometry expected with that data set. We simulate a realistic set of galaxy flux densities in the MIRI bands, using predictions from real galaxies observed in deep Hubble Space Telescope (HST) imaging from CANDELS (Grogin et al. 2011; Koekemoer et al. 2011; Stefanon et al. 2017). We then generate simulated raw MIRI images that account for realistic estimates of the noise. We then reduce MIRI raw imaging data to obtain MIRI mosaics that are matched to the existing HST/CANDELS imaging. We measure the mid-IR fluxes in the MIRI bands. We then perform SED fitting with x-cigale on the fluxes of MIRI and currently existing bands in UV to IR wavelengths. The redshift and other source properties are fit simultaneously in this SED-fitting process, where results are then derived by marginalizing over the other parameters. Finally, we evaluate the results by comparing the SED-fitting source properties with the model input ones.

This paper is structured as follows. In Section 2, we describe the CEERS survey, the existing CANDELS catalog, the construction of our model MIRI flux densities, and how we generated the simulated MIRI (raw) data. In Section 3, we discuss the reduction of the simulated MIRI raw data, discuss the method for measuring matched-aperture flux densities, and analyze the SEDs (using SED fitting from x-cigale). We then assess the results, focusing on the ability to constrain the photometric redshifts and fracAGN of distant galaxies by including the MIRI data. In Section 4, we discuss our results and future prospects using JWST/MIRI for this application. We summarize our work in Section 5.

Throughout this paper, we assume the same cosmology as x-cigale, i.e., a flat ΛCDM cosmology with H0 = 70.4 km s−1 Mpc−1 and ΩM = 0.272 (WMAP7; Komatsu et al. 2011). We adopt a Chabrier initial mass function (IMF; Chabrier 2003) for relevant quantities (stellar masses, etc.). Quoted uncertainties are at the 1σ (68%) confidence level. All magnitudes are in AB units (Oke & Gunn 1983), where ${m}_{\mathrm{AB}}=-48.6-2.5\mathrm{log}({f}_{\nu })$ for fν in units of erg s−1 cm−2 Hz−1.

2. Sample and Data2.1. The CEERS Survey

The CEERS Survey is an approved JWST program (Finkelstein et al. 2017), covering ≈100 arcmin2 in the Extended Groth Strip field (EGS). CEERS consists of 63 hr of NIRCam (1–5 μm) and MIRI (5–21 μm) imaging, NIRSpec R ∼ 100 and R ∼ 1000 spectroscopy, and NIRCam/grism R ∼ 1500 spectroscopy. CEERS will be one of the first public JWST surveys with data publicly available 5 months after acquisition in Cycle 1.

CEERS has four MIRI pointings, labeled as MIRI 1–4 fields, respectively. In this paper, we focus on the MIRI2 pointing that has the widest wavelength coverage of six bands: F770W, F1000W, F1280W, F1500W, F1800W, and F2100W. MIRI2 is centered at (214fdg95330, 52fdg95101).16 For the five bands from F770W to F1800W, the planned exposure time is 1665 s per band (three dithers per band); for F2100W, the exposure time is 4662 s (six dithers). The designed 5σ limiting magnitudes in the proposal are 25.5 (F770W), 24.8 (F1000W), 24.3 (F1280W), 23.8 (F1500W), 22.9 (F1800W), and 22.8 (F2100W). Figure 1 displays the transmission curves of the available MIRI bands in the MIRI2 pointing.

Figure 1.

Figure 1. Top: star-forming galaxy (black) and AGN (red) SEDs at z = 1.5 generated by x-cigale. The prominent PAH emission features are marked by the vertical lines (at rest frame 6.2, 7.7, 8.6, and 11.2 μm; see Tielens 2008). Bottom: transmission curves of the MIRI filters for the CEERS/MIRI2 pointing (adopted in this work). By covering the PAH features and AGN hot-dust emission, the MIRI data can potentially improve the constraints on photo-z and the AGN component.

Download figure:

Standard image High-resolution image

With the goal of making our simulation as realistic as possible, we use all available data for the galaxies in the CEERS/MIRI2 field. Our work is based on the F160W (H-band) selected CANDELS/EGS catalog (Stefanon et al. 2017). In this catalog, there are 463 sources down to H160 ≈ 26.6 within the MIRI2 pointing. Only 8 of these sources have secure spectroscopic redshift (spec-z) measurements, while the rest have photometric redshifts (photo-z) from Stefanon et al. (2017). These CANDELS/EGS sources are mostly in the redshift range of z = 0–3. The CEERS/MIRI2 region is covered by 15 broad bands from CFHT/MegaCam u* to Spitzer/MIPS 24 μm (Dickinson et al. 2006; Stefanon et al. 2017). Although this region is also covered by Herschel (Lutz et al. 2011; Oliver et al. 2012), most of our sources are beyond the sensitivity of Herschel (only two sources in the current MIRI2 field are detected (>5σ) by PACS, and none by SPIRE). Therefore, we do not use Herschel data in this work owing to the low detection rate. Based on these photometric and redshift data, we obtain the model input MIRI fluxes in Section 2.2.

2.2. The Model MIRI Fluxes

We have incomplete knowledge as to which galaxies in the CANDELS catalog (within MIRI2 pointing) host AGNs, because the currently existing photometric data provide very weak constraints on AGN emission (see Section 3.3.2). However, considering the relatively small area of the MIRI FOV (2.2 arcmin2), we estimate that only ∼10 AGNs (out of 463 sources; Section 2.1) actually exist within the MIRI2 pointing, based on the AGN number density derived from the deepest X-ray survey, the 7 Ms Chandra Deep Field-South (CDF-S; Luo et al. 2017). Actually, none of the MIRI2 galaxies have been detected by the existing 800 ks X-ray data in the EGS field (2–7 keV; Nandra et al. 2015). Therefore, in Section 2.2.1, we first model the SEDs of real existing photometry with pure-galaxy models (i.e., fracAGN = 0) and integrate the best-fit SEDs with the MIRI filter transmission curves to obtain the bandpass-averaged flux densities as inputs to the MIRI simulations. To test the cases when an AGN is present, we also add a hypothetical AGN component to the best-fit galaxy SEDs and rederive the MIRI flux densities to test for this possibility in additional sets of simulations in Section 2.2.2. These MIRI fluxes are used as the model input in our simulations of MIRI imaging data in Section 2.3.

2.2.1. Pure-galaxy Models Based on Empirical SEDs

We obtain model MIRI fluxes based on the SED fitting of real observations. We fit the currently existing broadband photometric data (from Stefanon et al. 2017) that exist for the galaxies in the nominal MIRI2 pointing (Section 2.1) with x-cigale (Boquien et al. 2019; Yang et al. 2020). We do not use the X-ray module of x-cigale owing to the lack of X-ray detections (Section 2.2). In the x-cigale analysis, we fix the redshift at the spec-z value (when available) or otherwise the photo-z from the CANDELS/EGS catalog (Stefanon et al. 2017). The fitting parameters are similar to those used by Yang et al. (2020) and are summarized in Table 1. As in Table 1, we adopt the model of Dale et al. (2014) for galactic dust emission. In this model, a single parameter of radiation α slope17 controls the IR SED shape (e.g., the MIR/FIR ratio and PAH emission strength). The model is derived from observations of local galaxies, and admittedly it may deviate from the true IR SEDs of distant sources. However, for this exercise it is sufficient, as our goal is to recover the mid-IR emission. In the future, new models dedicated for distant galaxies can be obtained using empirical measurements from the MIRI data themselves (especially if additional spectroscopy from the medium-resolution spectrometer, MRS, becomes available). We will implement these new models to x-cigale to facilitate the modeling of distant galaxies with real MIRI data.

Table 1.  x-cigale Parameters for the Fitting of Existing Data

Module Parameter Values Star formation history τ (Gyr) 0.1, 0.5, 1, 5 $\mathrm{SFR}\propto t\exp (-t/\tau )$ t (Gyr) 0.5, 1, 3, 5, 7 Simple stellar population IMF Chabrier (2003) Bruzual & Charlot (2003)     Stellar extinction     Calzetti et al. (2000) E(B − V) 0.1, 0.2, 0.3, 0.4, Leitherer et al. (2002)   0.5, 0.7, 0.9 Galactic dust emission Radiation 1.0, 1.5, 2.0, 2.5 Dale et al. (2014) α slope   Redshifting z Fixed at z of     Stefanon et al. (2017)

Note. For parameters not listed here, we use x-cigale default values.

Download table as:  ASCIITypeset image

After running x-cigale, we obtain the model MIRI fluxes by integrating the best-fit SEDs with the transmission curves of MIRI filters (Figure 1). Figure 2 shows the best-fit dust IR luminosity as a function of redshift. Because we do not have FIR data to directly constrain galaxy cold-dust emission, the IR luminosity largely comes from the modeling of UV/optical stellar extinction, but this is satisfactory for our purposes here, where we study the measured mid-IR emission from MIRI based on the simulations. x-cigale follows energy conservation, and thus the extincted UV/optical luminosity equals the dust reemitted IR luminosity.

Figure 2.

Figure 2. Total dust IR luminosity vs. redshift for our simulated sources in the CEERS/MIRI2 pointing. The redshift distribution is displayed in the top panel. The IR luminosities are obtained from x-cigale SED fitting with pure-galaxy templates (Section 2.2.1). The SFR is labeled on the right, which is converted from LIR using Kennicutt (1998) modified for our adopted Chabrier IMF (Salim et al. 2007). The sources detected (>5σ) in at least two MIRI bands are highlighted as blue squares. As expected, the MIRI-detected sources tend to be more IR luminous than the undetected ones. The orange circles mark the sources detected (>5σ) by MIPS 24 μm. The CEERS/MIRI data will probe galaxies with IR luminosities an order of magnitude lower than the current MIPS data (at fixed redshift).

Download figure:

Standard image High-resolution image 2.2.2. Hypothetical AGN–Galaxy Mixed Models

In Section 2.2.1, we derive the model input MIRI fluxes from pure-galaxy templates based on the real existing photometry. Now, we add a hypothetical AGN component to the SEDs to explore the cases when an AGN is present. The purpose of using realistic galaxy templates is to test, if a galaxy in our field hosts an AGN, how well we can detect/characterize the AGN component. As we argue below (Section 4.2), MIRI observations may be sensitive to the presence of emission from obscured AGNs even in galaxies that are not detected in the X-ray data. For this reason, it is useful to test the ability to recover the AGN emission in these "real" objects like those that will be seen in CEERS.

To generate AGN–galaxy mixed SEDs, we again employ x-cigale, which not only fits the photometric data but also generates model SEDs with given physical parameters. For each source, we adopt the best-fit parameters in Section 2.2.1 for the galaxy component. Then, we add an AGN component with the SKIRTOR model in x-cigale (Stalevski et al. 2012, 2016). SKIRTOR is a modern clumpy torus model based on Monte Carlo simulations. Yang et al. (2020) introduced a PD component to this model when implementing it to x-cigale.

For the AGN model parameters, we set the torus angle (between horizontal plane and torus edge) as 40° (the x-cigale default value), which is favored by the observations of Stalevski et al. (2016). We randomly choose a viewing angle (between AGN axis and line of sight) among 60°, 70°, 80°, and 90°. We do not include smaller face-on viewing angles that correspond to type 1 (broad-line) AGNs. These AGNs are typically luminous, as the UV/optical emission from the central engine is directly visible. Therefore, their properties such as redshift and AGN luminosity can be well measured with UV/optical spectroscopy and photometry. Also, most AGNs selected in deep X-ray surveys are type 2 (see Section 1). In this paper, we focus on type 2 AGNs for which the UV/optical emission is entirely obscured, and therefore our choice of the larger viewing angles is appropriate.

For other AGN torus parameters such as the 9.7 μm optical depth (τ9.7) and radial profile index p (density ∝ r−p), we randomly assign values to each source from all allowed values in the SKIRTOR model (e.g., τ9.7 = 3, 5, 7, 9, 11 and p = 0, 0.5, 1, 1.5). For the AGN PD E(B−V), temperature, and emissivity, we randomly assign values for each source from ranges of 0–0.3 (SMC extinction law), 100–300 K, and 1–2. This random choice of model parameters is to represent the diversity of torus physical properties based on current constraints (Stalevski et al. 2012, 2016). There are a total of nine parameters describing the hypothetical AGN component. We note that four of these parameters are free when we perform the SED fitting, while the other five parameters are fixed to canonical values (see Section 3.3). We further note that the values of these fixed parameters do not impact the conclusions here (see also Yang et al. 2019).

After determining the AGN parameters above, we generate the total (galaxy + AGN) SEDs with a given fracAGN and convolve the SEDs with MIRI filter curves to obtain model input MIRI fluxes. We adopt five fracAGN values of 0.2, 0.4, 0.6, 0.8, and 0.99, for different AGN strengths.18 For each fracAGN, we perform a set of MIRI simulations in Section 2.3. Figure 3 shows a set of SED models with different fracAGN values (as a reminder, fracAGN is the fraction of the IR luminosity [3–1000 μm] from the AGN component). As fracAGN increases, the shape of mid-IR SEDs changes significantly, because AGN hot-dust emission mostly concentrates on mid-IR wavelengths. After inserting the hypothetical AGN component, the flux densities at other bands may exceed the observed fluxes (this is especially the case for the IRAC and MIPS bands), and we deal with this issue in Section 3.3 by generating mock fluxes for these bands from the new SED.

Figure 3.

Figure 3. Set of x-cigale SED templates with different fracAGN. The Fν is normalized at 600 nm. The gray shaded region indicates the MIRI coverage (F770W to F2100W) at z = 1.5. The type 2 AGN component mainly affects the SED of mid-IR wavelengths, which are covered by MIRI. MIRI also covers the galactic PAH emission features.

Download figure:

Standard image High-resolution image 2.3. The Simulations of Imaging Data

We use mirisim (version 2.2.0; Klaassen et al. 2021) to simulate the images of MIRI bands. mirisim, developed by the MIRI European Consortium, is a dedicated simulation package for the JWST/MIRI instrument. mirisim accepts source positions, fluxes, and morphological shapes as input. The MIRI images generated by mirisim are preliminary (raw) "level 0" data products and need to be processed through the jwst calibration pipeline19 (pipeline hereafter) for further reduction (Section 3.1). In this way, mirisim provides data products equivalent in format to what we expect JWST to deliver.

We set the background noise as "low" in mirisim, as the EGS region has weak zodiacal background (Finkelstein et al. 2017). We set the background spatial gradient to 5% per arcminute (i.e., the background changes by 5% per arcminute) with a position angle (PA) of 45°. We adopt a dither pattern identical to the one in the CEERS Astronomer's Proposal Tool (APT) file, optimized for parallel observations between the Near Infrared Camera (NIRCam) and MIRI (Finkelstein et al. 2017). The exposure duration for each dither is the same as that in the CEERS APT configuration (using the same Groups/integration, integrations/exposure, and exposures/dither as expected for CEERS; Section 2.1). The detector readout mode is set to "fast," the same as proposed for CEERS. The input fluxes are the model fluxes obtained in Section 2.2.

Besides a point-like source profile, mirisim also allows galaxy morphologies with Sérsic profiles. Based on the HST F160W imaging data, van der Wel et al. (2012) performed Sérsic fitting (Sersic 1968) for H160 ≤ 24.5 galaxies in the CANDELS/EGS field. We adopt their results for sources with good fitting quality (flag ≤ 1). For the rest of the sources (flag > 1 or H160 > 24.5), we adopt a point-like profile as mirisim input to test results for unresolved sources. The adopted galaxy profiles have diameters (2 × Re) of 0farcs27–0farcs85 (20th–80th percentile), while the MIRI PSFs have FWHMs ranging from 0farcs24 (F770W) to 0farcs67 (F2100W). Therefore, MIRI can spatially resolve a nonnegligible fraction of the CANDELS galaxies even at ≳20 μm. In this work, we assume that the AGN and galaxy emissions have the same morphology, but in reality the former is likely more compact than the latter. Therefore, a source could be more compact in MIRI than in F160W owing to the presence of an AGN. Such a morphological difference between MIRI and F160W could serve as an AGN indicator. We leave the investigation of this potential technique of AGN identification to future works.

For each model input fracAGN (Section 2.2), we perform an independent set of mirisim simulations for all six MIRI bands (from F770W to F2100W). Figure 4 displays the simulated images of the pure-galaxy case in false (RGB) color (after the data reduction in Section 3.1). It is clear that the colors and morphologies from MIRI will provide an unprecedented view of the IR emission from distant galaxies compared to previous missions (e.g., MIPS 24 μm).

Figure 4.

Figure 4. Images of HST/WFC3 F160W, JWST/MIRI F770W+F1000W+F1280W (blue–green–red false color), JWST/MIRI F1500W+F1800W+F2100W (blue–green–red false color), and Spitzer/MIPS 24 μm. The F160W and MIPS images are real; the MIRI images are simulated. The four zoom-in boxes have a size of 10'' × 10'' and highlight various examples of galaxies with different MIRI colors and levels of crowding. The MIRI cutouts are from the simulation set of pure-galaxy models (Section 2.2.1). Different MIRI colors reveal different observed-frame wavelengths of the PAH features due to redshifting (e.g., Figures 1 and 3). Note that the angular resolution of MIRI is much higher than that of MIPS and enables detections of galaxies previously blended by MIPS.

Download figure:

Standard image High-resolution image

Since our simulations are based on the CANDELS F160W sources, we assume that all MIRI sources are detected in F160W. We acknowledge that a new type of optically faint mid-IR-bright sources could actually exist, and such sources may be detectable by MIRI but not by F160W. The real MIRI images after launch will provide a great chance to explore this unknown regime of distant galaxies.

In addition to the real CANDELS sources, we also simulated a set of bright point sources in Appendix A using the CEERS MIRI exposure configuration. These point sources are used for the purposes of data validation, photometric tests, and corrections (see Appendix A). For the photometry extraction purpose (see Section 3.2.1 for details), we also need to simulate HST F160W images using the identical galaxy morphologies (i.e., pure Sérsic models) that are equivalent to those in the simulated MIRI images. We describe the simulated F160W images in Appendix B.

3. Data Processing and Analyses3.1. The Reduction of MIRI Imaging Data

We run the JWST pipeline (version 0.15.2) to process the mirisim output of the level 0 data from mirisim (Section 2.3). The level 0 data for each dither are in the format of a FITS file. There are three stages of data processing. Stage 1 performs detector-level corrections for individual dithers, producing count-rate maps from the "up-the-ramp" readouts. Stage 2 applies physical corrections and calibrations to individual dithers, producing flux maps with image header containing standard World Coordinate System (WCS) information. Stage 3 combines individual dithers into a single image.

Stage 1 involves the rejection of cosmic rays, which are deposited onto the imaging data during the exposure. While this process successfully rejects most (>99%) of the cosmic rays, it still misses up to ≈20 cosmic rays per dither after visual inspection. To clean up these outliers, we run astro-scrappy (McCully et al. 2018), which implements the "Laplacian Edge" algorithm of cosmic-ray detection (see van Dokkum 2001), and mask the pixels polluted by the detected cosmic rays. Stage 3 involves background subtraction before merging the dithers. However, the current version of pipeline only subtracts a single global background value, although the actual background may have spatial dependence due to instrumental thermal emission and/or scattered light. Therefore, before merging the dithers, we perform an additional background subtraction with sep (version 1.0.3; Barbary 2016), which realizes the core algorithms of source extractor (Bertin & Arnouts 1996) in python. We adopt a cell size of 32 pixels (pixel size = 0farcs11) and a median-filtering size of 3 cells. Figure 5 demonstrates the effect of our background subtraction. After this background subtraction, we find that there is a small fraction of pixels (≲0.1%) that have extremely negative values (5σ F1500W sources). As we do not know what the morphologies of MIRI sources will be, we defer more detailed analysis to a future work.

Figure 8.

Figure 8. Left: tphot-measured F1500W AB magnitude vs. the input model magnitude in the mirisim simulation (Section 2.3). The black dashed line represents the equality relation (measured = input). Middle: measured minus input F1500W magnitudes (Δm1500) vs. the input magnitude. The red curve represents the running median of the data points with 20 sources per bin. The black dashed line indicates an offset of zero. Right: absolute value of σ1500 = ∣Δm1500∣ vs. the input magnitude. The red curve represents our error function. The black dashed vertical line represents our estimated 5σ limiting magnitude, where f/Δf = 5 from our error function. F1500W is shown here as an example, and the full version of this plot including all MIRI bands can be found in Appendix D.

Download figure:

Standard image High-resolution image Figure 9.

Figure 9. Difference between measured and input F1500W magnitudes (Δm1500) vs. the input magnitude, separated by Sérsic index. The red curves represent the running median. The orange circles highlight the sources (m1500,src) that have at least one neighbor ("nb," within 2'' radius) brighter than m1500,src + 1 mag. The top panel is for Sérsic n < 2 and point sources, and the bottom panel is for n > 2 sources. For the former, the measured magnitudes do not systematically deviate from the input magnitudes, but this is not true for the latter. This systematics for n > 2 sources is likely due to their extended faint wings. The full version of this plot including all MIRI bands can be found in Appendix D.

Download figure:

Standard image High-resolution image

In Figure 9, we indicate crowded sources (m1500,src) having at least one neighbor (within 2'' radius) brighter than m1500,src + 1 mag. These are objects where flux residuals from the (relatively) bright companion can bias the flux measurements of the sources. From Figure 9, the measured photometry of these sources is not always deviant from the distribution, although we notice that these sources do have a larger scatter (σMAD = 1.48 × MAD, where MAD is the median absolute deviation) in Δm1500 compared to isolated sources (0.23 vs. 0.17 for the >5σ sources in the top panel of Figure 9). The sources with neighbors only contribute to a small fraction of the >5σ population (e.g., 15% for F1500W). Therefore, we conclude that source crowding does not significantly affect the quality of our measured photometry, and it is clearly superior to previous instruments thanks to the high angular resolution of MIRI (see, e.g., Figure 4).

The photometric errors in the tphot output are based on the χ2 fitting of low-resolution source profiles (Section 3.2). We find that these errors are often significantly smaller (up to a factor of ≳10) than the actual differences between the measured and input magnitudes. This underestimation of uncertainties could be possibly due to the assumption that pixel noise is uncorrelated in tphot (Merlin et al. 2015, 2016). We do not adopt the tphot uncertainties. Instead, we establish an empirical "error function" for each MIRI band based on the distribution of input to measured flux densities. First, we bin our sources on their input MIRI magnitudes with each bin containing 30 sources. For each bin, we calculate the median input magnitude and the median difference between the measured and input magnitudes. Finally, we obtain the error function by linearly interpolating the median difference as a function of median magnitudes (see Figure 8, right). Although the error function for each band is derived based on the simulation set of pure-galaxy models (Section 2.2.1), we also apply it to the other simulation sets with nonzero fracAGN (Section 2.2.2). From the right panel of Figure 8, the uncertainty is ≈0.05 mag at the bright end and rises toward faint sources as expected. For each source, we estimate the magnitude uncertainty by evaluating the error function on the measured magnitude.

We estimate the 5σ limiting magnitude for each band as the value where the error function equals 0.217 mag, corresponding to δf/f = 0.2 (assuming standard error propagation, δm = 1.086 × δf/f). Figure 10 compares these measured limiting magnitudes and those expected for CEERS (Finkelstein et al. 2017) based on the JWST Exposure Time Calculator (ETC), assuming a point-source profile. From Figure 10, our measured 5σ depths are similar to those in the CEERS proposal (differences < 1 mag), where we interpret offsets as a result of limited sample sizes, our more realistic treatment of source surface brightness profiles and source crowding, and different background noise assumptions of mirisim (spatially dependent; Section 2.3) versus ETC (uniform). We consider the last factor as mainly responsible for the relative large difference between the F2100W depths of mirisim and ETC (0.95 mag; Figure 10), as the F2100W has the strongest background in our filter set. However, whether the background is spatially dependent or not can only be determined from the real MIRI imaging data after launch. For now, we caution that the ETC signal-to-noise ratio (S/N) could be overoptimistic for F2100W (and also F2550W, for which background is even stronger than for F2100W), because ETC does not account for the spatial dependence of background.

Figure 10.

Figure 10. The 5σ limiting magnitudes for different bands derived from our empirical error function. The dashed line represents limiting magnitudes for the CEERS configuration (Finkelstein et al. 2017) derived from the JWST/ETC. Our measured limiting magnitudes are similar to the expected limits, where we interpret offsets as being a limitation of the relatively small sample.

Download figure:

Standard image High-resolution image 3.3. The SED Fitting

We again use x-cigale to perform SED fitting of the galaxies' (simulated) photometry that includes measurement errors on the photometry. To assess the effects of MIRI data, we focus on the sources having significant (>5σ) detections in at least two MIRI bands (Figure 2). For the simulation set of pure-galaxy models (Section 2.2.1), these MIRI sources consist of 53% and 75% of the H160 < 26 and H160 < 25 population, respectively. For the simulation sets of AGN–galaxy models (Section 2.2.2), the detection rates are even higher owing to the additional AGN contribution. In contrast, the MIPS 24 μm detection rates are only 5% (H160 < 26) and 9% (H160 < 25) in the real data. The high detection rates of MIRI highlight its superior flux sensitivity.

For the currently existing non-MIRI bands, we cannot use the actually observed photometry directly in the SED fitting. The reason is that the observed photometry is not entirely consistent with that expected from the SED models used to generate model input MIRI fluxes (Section 2.2). This is especially a serious issue for the models of high fracAGN. When we add a strong hypothetical AGN component to the galaxy SED (Section 2.2.2), the model-expected IRAC and MIPS fluxes could be significantly higher than the observed values. Another issue is that our adopted model input redshifts are mostly photo-z from the CANDELS/EGS catalog (Section 2.2), and these photo-z are just an approximation, not the "true" values. However, to assess the accuracy of our SED-fitting results, we need the true redshifts instead of approximated ones (Sections 3.3.2 and 3.3.3).

To deal with the issues above, we use a mock catalog for each simulation set (Section 2.2), which is built following the technique in Boquien et al. (2019). First, we convolve the model SED with the filter transmissions to obtain the model flux for each source. We then perturb the model flux with a Gaussian fluctuation (σ = observed flux error in the CANDELS/EGS catalog) for each non-MIRI band. We adopt the perturbed fluxes in the x-cigale input, while keeping the original flux uncertainties. The mock catalog assumes that the model input SED parameters such as redshift and fracAGN are the "true" values solving the aforementioned issue, while it also perturbs the photometry to account for observational uncertainties. We note that our use of the mock, model photometry is only to gauge the ability of our methods to recover input parameters in this work. For real MIRI data after JWST launch, the observed non-MIRI photometry rather than the mock one should be used.

As noted previously, one benefit of x-cigale is that redshift and other source properties can be modeled simultaneously, although users can also fix the redshift at a given value such as spec-z. This feature is extremely valuable for deep-field sources, whose spec-z are challenging to measure.

The x-cigale run time for our simulated MIRI pointing (a few hundred sources) is ≲1 hr on a typical desktop/laptop. The run time will still be acceptable (≲1 day) even for a large sample of a million sources, thanks to the efficient parallel algorithm of x-cigale (Boquien et al. 2019).

3.3.1. x-cigale Configurations

The x-cigale parameters for our fittings are summarized in Table 2. The model parameters for the galaxy properties are the same as in Table 1. For the AGN component, we do not explore the full parameter space (nine parameters; Section 2.2.2), because many parameters (such as radial profile index and PD emissivity) only have limited effects on the AGN SED. Also, adopting all the possible parameters would take too much unnecessary computational time for practical purposes. Therefore, we only allow multiple values for four key parameters (i.e., fracAGN, τ9.7, PD E(B−V), and PD temperature) that significantly affect the AGN SED, while fixing other parameters at the single default value. Also, we allow the redshift to be a free parameter in x-cigale over a grid of z = 0.01–6 (in steps of ${\rm{\Delta }}\mathrm{ln}(1+z)=0.03$). Therefore, the Bayesian analysis of x-cigale also considers the uncertainties of redshift when analyzing other source properties such as fracAGN, because the analysis is performed on the entire parameter space. There are a total of 10 free parameters in our fitting (4 for galaxy, 4 for AGN, 1 for redshift, and 1 for SED normalization). We have 21 photometry data points (15 real existing bands and 6 simulated MIRI bands).

Table 2.  x-cigale Parameters for the Fitting of Simulated Data

Module Parameter Values Star formation history τ (Gyr) 0.1, 0.5, 1, 5 $\mathrm{SFR}\propto t\exp (-t/\tau )$ t (Gyr) 0.5, 1, 3, 5, 7 Simple stellar population IMF Chabrier (2003) Bruzual & Charlot (2003)     Stellar extinction     Calzetti et al. (2000) E(B−V) 0.1, 0.2, 0.3, 0.4, Leitherer et al. (2002)   0.5, 0.7, 0.9 Galactic dust emission Radiation α slope 1.0, 1.5, 2.0, 2.5 Dale et al. (2014)     AGN (UV-to-IR) Torus τ9.7 3, 5, 7, 9, 11   Torus angle 40° SKIRTOR Viewing angle 70°   fracAGN 0–0.99 (step 0.1)   PD extinction law SMC   PD E(B−V) 0–0.3 (step 0.1)   PD temperature (K) 100, 200, 300 Redshifting z 0.01–6 (step     ${\rm{\Delta }}\mathrm{ln}(1+z)=0.03$

Note. For parameters not listed here, we use x-cigale default values.

Download table as:  ASCIITypeset image

With the configurations above, we run x-cigale twice. First, we run x-cigale with the currently existing photometry only (Section 2.2). Second, we run x-cigale with the photometry from the mock catalogs (Section 3.3) and the simulated MIRI photometry (Section 3.2.2). We adopt the Bayesian (rather than best-fit) values of the source properties (such as redshift and fracAGN). Unlike the best-fit value, the Bayesian value properly considers all of the models weighted by their probabilities. The x-cigale runs are performed for the six simulation sets of pure-galaxy (Section 3.3.2) and galaxy–AGN mixed models (Section 3.3.3). Therefore, there are a total of 12 ( =2 × 6) x-cigale runs.

3.3.2. Fitting Results for Pure-galaxy Models

In this section, we present the SED-fitting results for the simulation set of realistic pure-galaxy models (i.e., fracAGN = 0; Section 2.2.1). Note that "pure galaxy" here refers to the input models that generate the MIRI photometry. In the fitting (Section 3.3.1), we allow both zero and positive fracAGN (Table 2), because we do not know whether an MIRI-observed galaxy hosts an AGN or not in realistic cases. We remind the readers that the output source properties (such as fracAGN) are continuously distributed, although the parameters in Table 2 are discrete. This is because the Bayesian results are pdf-weighted means of the parametric grid values (see Sections 3.3.1 and 4.3 of Boquien et al. 2019). As stated in Section 3.3, we focus on sources that have >5σ detections in ≥2 MIRI bands.

Figure 11 displays the fitting results for an example source. Note that the MIRI bands cover the PAH emission feature, which can be used as a robust redshift indicator and constrains the nature of the IR emission. Including the MIRI data greatly improves both redshift and fracAGN constraints, compared to the fitting results without MIRI. In particular, MIRI facilities a constraint on fracAGN, where data that lack MIRI provide almost no information. This is because MIRI covers the shape of the emission associated with the hot dust heated by the AGN (Section 1).

Figure 11.

Figure 11. Top: example SED fitting and residual of a source in the mock catalog. The galaxy has (true) input z = 2.13 and fracAGN = 0 (i.e., zero AGN component). The orange and purple curves indicate AGN and galaxy components, respectively. Since the best-fit fracAGN is 0 for this source, the galaxy component overlaps with the total SED. The MIRI data points are highlighted in red. The MIRI bands cover the PAH emission features, which are redshift indicators and constrain the nature of the IR emission. Bottom: 2D pdf and 1D pdf's of redshift and fracAGN for the source in the top panel. The 2D pdf density is scaled such that the integral over the 2D plane equals unity. The blue and red solid curves are from the fits with and without MIRI bands, respectively, as labeled. The dashed line indicates the Bayesian values from x-cigale output. The MIRI data are helpful in constraining redshift and fracAGN more tightly.

Download figure:

Standard image High-resolution image

In Figure 12, we compare the x-cigale redshift and fracAGN with the input values (Section 2.2). Following the convention in the literature (e.g., Yang et al. 2014), we adopt fractional redshift uncertainty, Δz/(1 + zinput), where Δz = zfit − zinput and zfit and zinput are x-cigale output and input model redshifts, respectively. We calculate the median and σMAD of Δz/(1 + zinput), as well as the outlier fraction (sources with ∣Δz∣/(1 + zinput) > 0.15). These results are marked in Figure 12. All three quantities (median, σMAD, and foutlier) are improved significantly after using MIRI. Notably, foutlier drops from 10% to 1.5% thanks to MIRI. One factor for the photo-z improvement is likely the capture of PAH emission by MIRI (e.g., Figure 1). We discuss the effectiveness of PAH emission in redshift constraint below.

Figure 12.

Figure 12. Top: measured source properties vs. the input source properties in the mirisim simulation, where we only use the galaxy component for the input MIRI fluxes (Section 2.3). The left and right panels are for redshift and fracAGN, respectively. The blue points and red squares are from the SED fits with and without MIRI photometry, respectively. In the left panel, the black dashed line indicates equality (input z = measured z), and the black dotted lines indicate a 15% redshift uncertainty. In the right panel, the x-axis positions are randomly perturbed for display purposes only, as the input fracAGN is 0. Bottom: distributions of (zfit − zinput)/(1 + zinput) (left) and fracAGN,fit − fracAGN,input (right). The median and σMAD values of these distributions are labeled. In the left panel, we also label the fraction of redshift outliers (∣Δz∣/(1 + z) > 0.15). With MIRI photometry, both of the uncertainties of redshift and fracAGN are smaller.

Download figure:

Standard image High-resolution image

First, we argue that the typical PAH emission is sufficiently strong to be detected by our MIRI filters. For example, the 6.2 μm line has a typical rest-frame equivalent width (EW) of ≈0.5 μm in our galaxy-dust models of Dale et al. (2014) (see, e.g., Armus et al. 2007; Spoon et al. 2007 for similar values). This EW translates to a flux excess of Δm1500 ≈ 0.4 for a z = 1.5 galaxy (Figure 1).20 This Δm1500 value is significantly larger than our >5σ sources' photometric uncertainties (≈0.05–0.22 mag; Section 3.2.2). The 7.7 μm PAH EW is typically ≈3 times wider than the 6.2 μm one (e.g., Stierwalt et al. 2014), and thus the former will even have a stronger impact on our MIRI photometry than the latter. We remind the reader that the observed EW grows with redshift as (1 + z), so the effects of PAHs on the MIRI bands are more substantive toward higher redshift. Therefore, the MIRI colors could provide useful information about PAH strength and photo-z. However, we acknowledge that a careful study is needed to quantify how useful the PAH features actually are in the JWST era. We will perform such a study in the future.

We also calculate the median and σMAD values of ΔfracAGN as shown in Figure 12. With MIRI, the median and σMAD of ΔfracAGN are 0.002 and 0.003, respectively, both close to 0. This result indicates tight constraints on the AGN component. Without MIRI, the constraints are poor, and the median and σMAD are 0.200 and 0.155, respectively, indicating that the constraints on the AGN component are poor. The significant role of MIRI in constraining AGNs highlights its ability of sampling the emission from AGN-heated hot dust (e.g., Figures 1 and 11). Without MIRI, there is a large gap in coverage from IRAC 8 μm to MIPS 24 μm, which is unable to probe the hot dust heated by the AGN, and thus the constraints on the AGN component are much weaker (e.g., Figure 11). Specifically, surveys that include multiple MIRI bands will be very effective at identifying galaxies without AGN emission.

Because the AGN hot-dust emission typically peaks at rest frame ∼10 μm (e.g., Figure 3), it is possible that the constraints on fracAGN become weaker at high redshifts, when the bulk of AGN emission shifts out of the MIRI coverage. To investigate this redshift dependence, we show the distributions of ΔfracAGN for different redshift bins in Figure 13. The constraints on fracAGN are excellent up to z ≈ 3, with median and σMAD both below ≈2%. At higher redshifts (z ≳ 3), the constraints become substantially weaker as expected, as the mid-IR features associated with hot dust from the AGN shift to higher wavelengths than can be probed by even MIRI. The photo-z quality also appears to drop at z ≳ 3, although the high-z sample size is not sufficiently large for a solid conclusion. The relatively poor photo-z quality at high z, if true, could result from the fact that the PAH features also shift to wavelength beyond those covered by MIRI.

Figure 13.

Figure 13. Same format as the bottom right panel of Figure 12, but dividing into different redshift bins as labeled. Values in blue ("w/ MIRI") correspond to results that include the MIRI data, and values in red ("w/o MIRI") exclude the MIRI data. With MIRI data, the constraint on fracAGN is robust at z ≲ 3. But the constraint becomes weaker at z ≳ 3, because the bulk of AGN hot-dust emission shifts out of MIRI coverage.

Download figure:

Standard image High-resolution image 3.3.3. Fitting Results for Galaxy–AGN Mixed Models

In this section, we present the SED-fitting results for the simulation sets of hypothetical galaxy–AGN mixed models with different positive values of input fracAGN (Section 2.2.2). Figure 14 displays an example SED fitting and pdf's for a source with model fracAGN = 0.4. After using MIRI data, the projected 1D pdf's of both redshift and fracAGN become narrower.

Figure 14.

Figure 14. Same format and source as Figure 11, but in the simulation with model input fracAGN = 0.4. Although the fracAGN Bayesian values (dashed lines) from fitting with and without MIRI are similar, the pdf of the fitting with MIRI is much narrower than that without MIRI.

Download figure:

Standard image High-resolution image

In Figure 15, we compare the distribution of redshift uncertainties from the fittings with and without MIRI data for different model input fracAGN. The redshift accuracy is improved after using MIRI data for each fracAGN case, as both the scatter and outlier fraction improve with the inclusion of MIRI data. The photo-z scatter, σMAD, increases as the input fracAGN increases. This is because, as AGN strength rises, the PAH features (used as a redshift indicator) become weaker in the SED (see Figure 3). Another reason is that we select more optically faint objects for higher fracAGN, as our selection is based on mid-IR (detected in ≥2 MIRI bands; Section 3.3). The optical flux uncertainties are larger for these optically faint sources. Regardless, there is appreciable gain when including the MIRI bands.

Figure 15.

Figure 15. Same format as the bottom left panel of Figure 12, but for the simulations with different model input fracAGN as labeled. Note that model fracAGN = 0.00 corresponds to pure-galaxy models (i.e., Figure 12 bottom). Values in blue ("w/MIRI") correspond to results that include the MIRI data, and values in red ("w/o MIRI") exclude the MIRI data. The redshift constraints are significantly tighter after adding MIRI data to the SED fitting.

Download figure:

Standard image High-resolution image

Figure 16 shows the fitting results of fracAGN. With MIRI, the dispersion of ΔfracAGN is remarkably small (≲0.003) for the cases of pure galaxy (fracAGN = 0) and AGN dominant (fracAGN = 0.99), and it is much larger (≳0.1) for the intermediate cases. This is expected, as in the two extreme cases (fracAGN = 0 or 0.99) the SED features of a galaxy or AGN are dominant. In contrast, when fracAGN is intermediate, the total SED has mixed features of a galaxy and AGN, and the SED decomposition is challenging. These sources would likely be considered "composites" in previous studies (e.g., Kirkpatrick et al. 2017).

Figure 16.

Figure 16. Distributions of fracAGN fitting results. From bottom to top, the model input fracAGN values are 0, 0.2, 0.4, 0.6, 0.8, and 0.99, respectively, as marked by the vertical black dashed lines. Values in blue ("w/MIRI") correspond to results that include the MIRI data, and values in red ("w/o MIRI") exclude the MIRI data. The fracAGN constraints become generally tighter after adding MIRI data to the SED fitting.

Download figure:

Standard image High-resolution image

At first glance, for model fracAGN = 0.4 and fracAGN = 0.6, it may appear that the accuracy of fracAGN does not change significantly after using MIRI data. However, this is misleading. This is because the Bayesian output value in x-cigale is calculated as the pdf-weighted mean (see Section 4.3 of Boquien et al. 2019), and this produces a median that tends to be located near the center of the parametric range (0.5 for fracAGN) in the case that the constraints are weak. Considering the extreme case when the model fracAGN = 0.5 and the pdf is totally flat, the Bayesian fracAGN will be exactly the same as the model value, although there are no constraints on fracAGN. The Bayesian fracAGN values (dashed lines) are similar for the fitting with and without MIRI, but the pdf's for individual objects are much narrower after using MIRI data. Therefore, sometimes it is not sufficient to use the mean values from the pdf only, and the errors (pdf-weighted standard deviation; Section 4.3 of Boquien et al. 2019) calculated by x-cigale may serve as a necessary diagnostic. Figure 17 displays the error distributions for all model fracAGN configurations. Indeed, for all model fracAGN (including fracAGN = 0.4 and fracAGN = 0.6), the median of fracAGN uncertainties always becomes smaller after using the MIRI data, showing that the constraint on fracAGN becomes tighter with MIRI. Quantitatively, MIRI data improve the fracAGN accuracy by a factor of ≈2 for the case of AGN–galaxy composite input (see Figure 17).

Figure 17.

Figure 17. Distribution of fracAGN errors from x-cigale. Different panels are for different model fracAGN as labeled. Values in blue ("w/MIRI") correspond to results that include the MIRI data, and values in red ("w/o MIRI") exclude the MIRI data. The blue and red colors indicate the SED fittings with and without MIRI data, respectively. The median values of the distributions are labeled in each panel. For each model fracAGN, the uncertainties become smaller after adding MIRI data to the SED fitting.

Download figure:

Standard image High-resolution image 3.3.4. Constraints on AGN Accretion Power with MIRI

In Section 3.3.3, we assess the SED-fitting quality of fracAGN for the cases where a hypothetical AGN is present in the input models. fracAGN describes the relative luminosities of an AGN versus a galaxy in terms of total IR luminosity. It is understandable that the fitted fracAGN still has significant scatter even using MIRI data, because we do not have far-IR photometry to tightly constrain the galaxy cold-dust emission, which also affects fracAGN.

However, it is often useful to obtain absolute AGN luminosities in AGN studies, as black hole (BH) accretion rates can be estimated from the absolute luminosities (e.g., Yang et al. 2019; Ni et al. 2019, 2021). x-cigale is able to estimate the intrinsic AGN accretion disk luminosity, Ldisk (i.e., "agn.accretion_power" in the output; Yang et al. 2020). The SKIRTOR AGN model in x-cigale adopts anisotropic disk emission, and Ldisk is calculated averaging over all viewing angles. Note that the Ldisk is numerically equivalent to the angle-averaged obscured disk + dust luminosity, as SKIRTOR is a physical model obeying energy conservation.

We compare the measured Ldisk and model input Ldisk in Figure 18 for the simulation sets of different input fracAGN (Section 2.2.2). The median and σMAD for all fracAGN cases are ≈−0.1 dex and ≈0.3 dex, respectively (see Figure 18). Therefore, we can reliably recover Ldisk from SED fitting of the photometric data of MIRI (and other bands). In the future, MIRI photometric surveys can be widely used in the studies of BH accretion and evolution across cosmic history (Section 4.2). From Figure 18, the quality of Ldisk becomes slightly better toward higher fracAGN, indicating that the constraint on AGN power is better for sources whose AGN SED component is more dominant in infrared.

Figure 18.

Figure 18. Measured vs. input AGN intrinsic disk luminosities. The measured Ldisk is based on the SED fitting with the simulated MIRI data, and different panels are for simulation sets of different input fracAGN as labeled. In each panel, the black dashed line represents the 1:1 relation, and the dotted lines represent ±0.5 dex from the 1:1 relation. The median and σMAD values of ΔLdisk (measured − input) are marked in each panel.

Download figure:

Standard image High-resolution image

Therefore, in addition to the constraints on fracAGN, the x-cigale results also provide reliable constraints on the BH accretion power, which is a useful quantity for AGN studies.

4. Discussion4.1. Comparison to Kirkpatrick et al. (2017)

Kirkpatrick et al. (2017) were one of the first to evaluate the use of MIRI to classify and characterize AGNs and star formation in distant galaxies. Their classification scheme is based on MIRI color–color diagrams and demonstrated that MIRI can identify objects whose IR emission is powered by star formation, AGNs, and composite cases. The method here uses x-cigale and SED fitting and improves on this previous work, as it uses all the available MIRI bands and measures simultaneously the photometric redshift and galaxy properties simultaneously, providing the pdf for each parameter.

We provide here a qualitative comparison between the two methods. Kirkpatrick et al. (2017) derived a parameter ${\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}$ (fractional AGN contribution to the rest-frame 5–15 μm SED, where the superscript K denotes that these are derived by Kirkpatrick et al. 2017). They then divided their source SEDs into three classes, i.e., galaxy (${\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}\lt 0.3$), composite ($0.3\leqslant {\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}\lt 0.7$), and AGN (${\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}\geqslant 0.7$), and define an empirical color–color scheme to classify simulated sources into these classes. They find that their color–color scheme is able to reach a high accuracy level of ≈79% (galaxy), 76% (composite), and 87% (AGN).

The thresholds ${\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}=0.3$ and ${\mathrm{frac}}_{\mathrm{MIR},\mathrm{AGN}}^{K}=0.7$ above roughly correspond to our fracAGN ≈ 0.1 and fracAGN ≈ 0.4, respectively (the exact conversions vary from model to model). Therefore, to compare with Kirkpatrick et al. (2017), we choose fracAGN < 0.1, 0.1 ≤ fracAGN < 0.4, and fracAGN > 0.4 as the criteria for the galaxy, composite, and AGN classes, respectively, in our analysis. Under these criteria, for our simulation set of pure-galaxy input models (Section 2.2.1), 88% of sources are correctly classified as a (star-formation-dominated) "galaxy" (i.e., a fracAGN in x-cigale output below 0.1; Figure 12). For the input models of fracAGN = 0.2 (i.e., the composite models; Section 2.2.2), the classification accuracy is 72% (Figure 15). For the input models of fracAGN = 0.6, 0.8, and 0.99 (i.e., the AGN models), the success rates are 77%, 90%, and 99%, respectively. Therefore, there is a high degree of overlap between the method of Kirkpatrick et al. (2017) and our SED-fitting method here.

Again, we note that the classifications above are for the comparison with Kirkpatrick et al. (2017) only. The method of Kirkpatrick et al. (2017) is straightforward: one just needs to apply a suitable color–color scheme depending on the source's redshift, which can be spec-z (if available) or photo-z. The qualitative classification results can be obtained instantaneously. Our SED fitting of x-cigale is a quantitative method that yields pdf-based property estimation and best-fit SEDs.

Another advantage of our x-cigale method is that the redshift does need to be known a priori, and x-cigale will fit simultaneously for the redshift and other parameters (including fracAGN). This feature is extremely valuable for deep-field sources, whose spec-z are challenging to measure.

The x-cigale run time for our simulated MIRI pointing (a few hundred sources) is ≲1 hr on a typical desktop/laptop. The run time will still be acceptable (≲1 day) even for a large sample of a million sources, thanks to the efficient parallel algorithm of x-cigale (Boquien et al. 2019).

4.2. Comparison with X-Ray AGN Selection

X-ray observations are effective in AGN identification (e.g., Brandt & Alexander 2015; Xue 2017). Strong X-ray emission is almost a universal property of the BH accretion process, and galactic processes (e.g., X-ray binaries and hot gas) can only reach low X-ray luminosities typically below LX ∼ 1042 erg s−1. Thanks to these strengths, X-ray surveys have detected numerous AGNs in the distant universe, significantly deepening our understanding of BH evolution across cosmic history (e.g., Civano et al. 2016; Luo et al. 2017; Yang et al. 2018a, 2018b). Like X-ray observations, MIRI can also reliably constrain the AGN accretion power (see Section 3.3.4). To evaluate the effectiveness of MIRI AGN selection, we compare the sensitivities of MIRI versus X-ray observations below.

To estimate the equivalent X-ray flux limit in the MIRI selections, first, we calculate model input AGN 6 μm luminosities (L6 μm) for all hypothetical sources in all of the fracAGN > 0 simulations (Section 2.2.2). We then convert L6 μm to LX based on the empirical L6 μm–LX relation in Stern (2015). We derive the X-ray fluxes (fX, observed-frame 0.5–7 keV) from LX assuming an X-ray spectral photon index of Γ = 1.8 (e.g., Yang et al. 2016; Liu et al. 2017). Finally, we decrease these fX by a factor of 2, which represents the typical X-ray obscuration effect (e.g., Luo et al. 2017). The fX distributions of MIRI-detected and undetected sources are displayed in Figure 19. From Figure 19, the MIRI detection becomes significantly incomplete (5σ) detected in at least two MIRI bands (for the default case where SEDs assume pure-galaxy models).3.  We perform x-cigale SED fitting, with and without the addition of the MIRI data (Section 3.3). We focus on the ability of the SED fitting to constrain the photometric redshift and fracAGN (these are fit simultaneously, and we marginalize the pdf's to derive constraints on these parameters). The accuracies of both redshift and fracAGN are improved significantly after using MIRI data, thanks to the capture of PAH features and AGN hot-dust emission by MIRI. Notably, for the simulation set of pure-galaxy models (which is likely the case for most MIRI sources), fracAGN can be constrained to the level of ≈0.2% with MIRI data, which is ≈100 times better than the fitting without MIRI data. At the same time, the photo-z scatter and outlier fraction are improved by a factor of ≈2 and ≈7, respectively, thanks to MIRI's capture of PAH emission.4.  Our SED-fitting method can reliably recover the source types (galaxy/composite/AGN; Section 4.1), similar to results from previous studies that were based on MIRI color–color diagnostics to characterize sources (e.g., Kirkpatrick et al. 2017).5.  We assess the AGN-detection sensitivity of MIRI versus X-ray (Section 4.2) and find that, for deep MIRI exposures (such as the CEERS-depth, 3.6 hr data simulated in this work), MIRI can already reach a sensitivity level even slightly higher than the deepest X-ray survey, CDF-S. MIRI should also detect Compton-thick AGNs for which X-ray selection is highly incomplete. Therefore, we conclude that MIRI will provide a complete census of the entire population of accreting BHs, enabling unbiased studies of BH evolution across cosmic history.6.  We discuss the effects of using different MIRI bands in SED fitting (Section 4.3). In particular, we focus on comparing the bluer (F770W, F1000W, and F1280W) and redder (F1500W, F1800W, and F2100W) bands. We find that the blue bands are more helpful for photo-z improvement than the red bands, because the former are more sensitive than the latter. However, the blue and red bands have similar effects in terms of constraining fracAGN, likely due to the fact that AGN emission is typically stronger in the red bands than in the blue ones.

Although this work focuses on the CEERS/MIRI2 observational strategy, it has general implications for other MIRI extragalactic surveys. For the surveys with the same filters (F770W to F2100W) but different depths, the qualitative conclusions in this paper are likely to hold in general, as any MIRI data with similar depth will achieve similar results. The only change is that different exposures will yield different detection limits (e.g., Figure 2). For example, deeper exposures will be able to detect more faint mid-IR sources at higher fidelity, thereby constraining their source properties (such as photo-z and fracAGN). For surveys with similar exposures to CEERS/MIRI2 but different band coverage, we have performed an example assessment in Section 4.3. We are willing to repeat the assessment for any other specific band set upon e-mail request, to satisfy the realistic needs of, e.g., survey planning.

We thank the referee for helpful feedback that improved this work. We thank the helpful discussions with Jacqueline Antwi-Danso, Emiliano Merlin, Karl Gordon, and the STScI JWST team. We also thank our collaborators on CEERS for their work and contributions to the project. The authors acknowledge the Texas A&M University Brazos HPC cluster and Texas A&M High Performance Research Computing Resources (HPRC, http://hprc.tamu.edu) that contributed to the research reported here. We also acknowledge our collaborators within CEERS for their input in the project. This work acknowledges support from the NASA/ESA/CSA James Webb Space Telescope through the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-03127. Support for program No. JWST-ERS-01345 was provided through a grant from the STScI under NASA contract NAS5-03127. P.G.P.-G. acknowledges support from Spanish Government research grant PGC2810-093499-BI00.

Software: astropy (v4.0.1; Astropy Collaboration et al. 2018), astro-scrappy (McCully et al. 2018), jwst calibration pipeline (v0.15.2; https://github.com/spacetelescope/jwst), mirisim(v2.2.0; Klaassen et al. 2021), sep (v1.0.3; Barbary 2016), source extractor (v2.19.5; Bertin & Arnouts 1996), tphot (v2.1; Merlin et al. 2015, 2016), x-cigale (Boquien et al. 2019; Yang et al. 2020).

Appendix A: Simulated MIRI Data Validation and Aperture Corrections

Both mirisim and pipeline are still in development. This work uses current versions of the available software, but we acknowledge that these could be updated prior to the acquisition of JWST data post-launch. One known issue is that mirisim and pipeline adopt different MIRI calibration files, which might also affect our simulated MIRI imaging data. Therefore, to investigate potential issues related to mirisim and pipeline, we applied additional tests to validate the imaging data before the photometry extraction process (Section 3.2). Below we perform data validation and correction processes based on a set of simulated bright point sources.

Using mirisim (Section 2.3), we simulate a grid of bright (but nonsaturated) point sources with an S/N of ≈1000 for each MIRI band. The separation between the neighboring sources is 15'', which is sufficiently large to avoid light pollution from neighbors. Figure 20 displays example simulated images for the MIRI F770W and F1800W bands (where these images have been fully reduced following the procedures in Section 3.1), where we have repeated this process for all MIRI bands to be obtained for the CEERS/MIRI2 field. These sources have a high S/N of 1000 and thus are less affected by the details of noise in the data reduction process compared to the realistic faint sources.

Figure 20.

Figure 20. Simulated point sources (after pipeline reduction) for the purposes of data validation and corrections. The left and right panels are for F770W and F1800W, respectively. The sources at different positions appear to have similar brightness, suggesting that the flat-field correction in pipeline is reliable.

Download figure:

Standard image High-resolution image

We measure the EEF as a function of radius for each simulated point source at different positions. The EEF at a given radius is calculated as the encircled flux within the radius divided by the model input flux. The results for F770W and F1800W are displayed in Figure 21 as examples. The EEFs for the point sources at different positions show slight variations, which is expected given the spatial variations in the background (Section 3.1), but the scatter is small (≈1% at 1''). This small scatter indicates that the flat-field correction in pipeline functions well. However, we find systematic offsets between the measured EEFs for these point sources compared to the theoretical PSF EEF from webbpsf (version 0.90).22 The specific reason is possibly due to the inconsistencies between the MIRI calibration files adopted by mirisim and pipeline (K. Gordon 2020, private communication). It is also possible that one of the simulation steps is truncating the wings of the PSF. As these possibilities are associated with the production of the simulations, they are unlikely to be present in the real JWST data. Here, for our purposes we measure the systematic offset as a "scaling correction," which we then apply to the flux densities in MIRI that we measure. This scaling "correction factor" is calculated as

Equation (A.1)

where med(EEFPT,0.8) is the median of the point-source EEF at r0.8 (the radius where PSF EEF = 0.8; see Figure 21). The value of 0.8 is chosen, because most of the light is encircled within r0.8 and the pixels are still source dominated instead of background dominated at r0.8. The values of fscaling are all within ≃10% of unity and range from 0.87 (F1800W) to 1.08 (F770W). We apply the scaling correction by multiplying all pixel values by a factor of fscaling for each MIRI band. After this procedure, the point-source EEFs become similar to the PSF EEF as displayed in Figure 21.

Figure 21.

Figure 21. EEFs for the simulated bright point sources of Figure 20 in F770W (top) and F1800W (bottom). The black and blue curves represent the EEFs before and after the scaling correction (Appendix A). The red solid curve is from the PSF generated by webbpsf. The red dashed vertical line indicates the radius corresponding to PSF EEF = 0.8, where the scaling correction factor is estimated.

Download figure:

Standard image High-resolution image

Besides the scaling correction above, we also need to consider aperture corrections. This is because the PSF faint wings extend to large radii and can contribute to the background (Figure 21). We calculate the aperture correction by comparing the total (input/true) flux density to the averaged measured flux density. The aperture correction factor is

Equation (A.2)

where EEFPSF,1farcs5 is the PSF EEF at 1farcs5. We choose the radius of 1farcs5 as a fiducial value, as beyond this radius the EEF has little growth (see Figure 21). We apply the aperture correction by multiplying the measured source fluxes by a factor of faper. As expected, faper is higher for longer wavelengths, which correspond to more extended PSF. The values of faper increase with PSF FWHM (which scales with wavelength) and range from 1.10 (F770W) to 1.19 (F2100W).

We note that the correction procedures above are not limited to the simulated MIRI data and can be applied to real MIRI data. After the launch of JWST, some bright "standard stars" with known mid-IR fluxes can play the same role as the simulated point sources in this work. In fact, some key calibration programs focusing on standard stars have already been scheduled for the upcoming Cycle 1.23

After the corrections above, we extract photometry for the simulated point sources using tphot (see Section 3.2). We assess the photometry quality by comparing the measured magnitudes and the model input magnitudes. Figure 22 shows the distribution of the magnitude offsets for F770W and F1800W as examples. The systematics and scatters are both ≲0.02 mag for all MIRI bands, indicating that our photometric corrections above and tphot measurements are reliable.

Figure 22.

Figure 22. Stacked-histogram distribution of Δmag (tphot-measured − model input) of F770W (top) and F1800W (bottom) for the simulated bright point sources of Figure 20 (where "fake pt" stands for "fake point sources"). The median and σMAD values are labeled. The small systematics and scatters indicate that our photometric corrections and tphot measurements are reliable.

Download figure:

Standard image High-resolution image Appendix B: The Simulation of F160W Imaging Data

To extract the photometry on the simulated MIRI images, we also need to simulate F160W images that have source morphologies consistent with MIRI morphologies (see Section 3.2.1). Note that we use the simulated HST image only to measure MIRI flux densities from TPHOT (we use the real F160W source fluxes and errors from the CANDELS/EGS catalog when fitting the SEDs in Section 2.2.1). First, we create an empty image covering the same region as the simulated MIRI fields, with pixel size set to 0farcs06 (the same as the CANDELS survey). We then position all sources (Stefanon et al. 2017) within the FOV on the image. The fluxes are set to the observed F160W fluxes in the Stefanon et al. (2017) catalog, and the morphologies use the same Sérsic model parameters (from the van der Wel et al. 2012 catalog) that are used for the mirisim input above. In addition, for point-like sources, we allow slight shifts to center the pixel at the pixel nearest to the source position; for a Sérsic source, we generate its profile utilizing the sersic2d function of astropy (Astropy Collaboration et al. 2018), which allows fractional pixel positions. Next, we convolve the image with the F160W PSF derived by the 3DHST team (Momcheva et al. 2016). We perform this convolution process with the convolve_fft function of astropy. Finally, we add a random Gaussian noise to each pixel with an amplitude from the rms map produced by the CANDELS team (Grogin et al. 2011; Koekemoer et al. 2011). Figure 23 compares the cutouts of the original/simulated F160W images and two MIRI-band images. For the galaxy near the center, the simulated and original F160W profiles are obviously different, since a Sérsic profile, as a proxy, cannot reproduce the complex morphological feature of spiral arms. As expected, the MIRI images have similar profiles to the simulated F160W image, because they share the same input Sérsic profile.

Figure 23.

Figure 23. Top: real (left) and simulated (right) HST F160W 30'' × 30'' cutouts with the same arcsinh color scale. The galaxy morphologies in the real and simulated images are broadly similar but not exactly the same. For example, the simulation does not reproduce the spiral arms of the galaxy at the center, since the simulation is limited to Sérsic profiles. Bottom: simulated (pure-galaxy models) MIRI F770W (left) and F1500W (right) images of the same sky region as the top panels. The MIRI morphologies are more consistent with that of the simulated F160W than the real F160W by design. This morphological consistency is the reason why we need to use the simulated F160W in photometry extraction.

Download figure:

Standard image High-resolution image Appendix C: Comparison between Source Extractor and Tphot Photometry

source extractor (Bertin & Arnouts 1996) is widely used for photometry extraction (e.g., Caputi 2013; Yang et al. 2014; Kauffmann et al. 2020). Therefore, it is instructive to test the performance of source extractor on our simulated MIRI images. Using F1500W as an example, we perform source extractor photometry and compare the resulting quality with that from tphot below. The results are similar for other MIRI bands.

The "AUTO" mode photometry in source extractor provides an estimation of the bulk of the total flux based on an adaptive-aperture algorithm (Kron 1980). The AUTO method automatically adjusts the size and shape of an elliptical aperture, and it works best for bright resolved sources with high S/Ns. However, for a faint low-S/N source, the aperture size (and thus the resulting flux) could be systematically underestimated (Bertin & Arnouts 1996). This issue can be alleviated by adopting a circular photometry aperture of a fixed size for faint sources. Therefore, we adopt a strategy of choosing the brighter one of the AUTO and fixed-aperture fluxes (after corrections) for each source.

We run source extractor (v2.19.5) on the same F1500W image as in our tphot run (Section 3.2.2; pure-galaxy inputs). Our source extractor run adopts the parameters listed in Table 3. These adopted parameter values are chosen to balance between detection purity and completeness. Using these parameters, we detect most (87%) of the input sources that are brighter than the ETC-estimated F1500W 5σ depth. We visually check the 13% undetected sources and find that most of them are either at the detector edges (where noise is strong) or blended with nearby bright sources. On the other hand, our adopted parameters result in a low false-detection rate (8%). Lowering the detection threshold would identify more of the missed sources (improving completeness) but also increase the number of "false detections" (which we can measure from our simulations, as we have the true, input catalog). Therefore, we consider that our source extractor configurations are balanced to be optimal. Our adopted aperture size (PHOTO_APERTURES = 31) corresponds to 80% point-source EEF in F1500W, and thus we correct the resulting aperture fluxes by a factor of 1/0.8 to estimate the total fluxes. Based on our choice of "PHOTO_AUTOPARAMS," the AUTO flux is expected to recover ≈94% of the total flux on average (Bertin & Arnouts 1996). Therefore, we correct the resulting AUTO fluxes by a factor of 1/0.94 to obtain the total.

Table 3.  Adopted source extractor Parameter Values for the F1500W Photometry

Parameter Value DETECT_MINAREA 9 DETECT_THRESH 2.3 FILTER Y FILTER_NAME Gauss_3.0_7x7.conv DEBLEND_NTHRESH 64 DEBLEND_MINCONT 0.00001 CLEAN Y CLEAN_PARAM 0.7 MASK_TYPE Correct PHOTO_APERTURES 31 PHOTO_AUTOPARAMS 2.5, 3.5 BACK_SIZE 64 BACK_FILTERSIZE 6 BACKPHOTO_TYPE LOCAL BACKPHOTO_THICK 40

Note. For parameters not listed here, we use default values.

Download table as:  ASCIITypeset image

Figure 24 compares the resulting source extractor magnitudes with the input magnitudes. Figure 25 compares source extractor versus tphot magnitudes for point/n < 2 and n ≥ 2 sources (similar to Figure 9). The systematic differences (as indicated by the median) between source extractor and tphot magnitudes are generally small (≲0.05 mag). This systematic consistency indicates that, like tphot, source extractor also misses the light from the faint extended wings for n ≥ 2 sources (see Section 3.2.2).

Figure 24.

Figure 24. Similar format to Figure 8, but for source extractor photometry. In the right panel, we also plot the tphot error function for comparison. The source extractor uncertainties are larger than tphot uncertainties in general.

Download figure:

Standard image High-resolution image Figure 25.

Figure 25. Difference between source extractor and TPHOT F1500W magnitudes (Δm1500) vs. the source extractor magnitude. The red curves represent the running median. The top panel is for Sérsic n < 2 and point sources, and the bottom panel is for n > 2 sources. For n ≥ 2 sources, the median Δm1500 is ≈ 0, indicating that, like tphot, source extractor also misses the light from the faint extended wings (see Section 3.2.2).

Download figure:

Standard image High-resolution image Appendix D: Complete Comparisons between Input and Measured Flux Densities for All MIRI Bands

In this appendix, Figure 26 shows the comparison between input and measured flux densities for all the MIRI bands (and is similar to Figure 8 above). Figure 27 compares the results dividing the sample by morphology (Sérsic index) for all bands (and is similar to Figure 9 above).

Figure 26.

Download figure:

Standard image High-resolution image Figure 26.

Figure 26. Same format as Figure 8, but for all MIRI bands.

Download figure:

Standard image High-resolution image Figure 27.

Figure 27. Same format as Figure 9, but for all MIRI bands.

Download figure:

Standard image High-resolution image


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3